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Abstract. We consider a two-scale reaction diffusion system able to capture the corrosion 
of concrete with sulfates. Our aim here is to define and compute two macroscopic corrosion 
indicators: typical pH drop and gypsum profiles. Mathematically, the system is coupled, 
endowed with micro-macro transmission conditions, and posed on two different spatially- 
separated scales: one microscopic (pore scale) and one macroscopic (sewer pipe scale). We 
use a logarithmic expression to compute values of pH from the volume averaged concentration 
of sulfuric acid which is obtained by resolving numerically the two-scale system (microscopic 
equations with direct feedback with the macroscopic diffusion of one of the reactants). Further- 
more, we also evaluate the content of the main sulfatation reaction (corrosion) product — the 
gypsum — and point out numerically a persistent kink in gypsum's concentration profile. Fi- 
nally, we illustrate numerically the position of the free boundary separating corroded from 
not-yet-corroded regions. 



1. Introduction 

1.1. Background on sulfate corrosion. Often in service-life predictions of concrete structures 
(e.g., sewer systems), the effects of chemical and biological corrosion processes are fairly neglected. 
In sewer systems and wastewater treatment facilities, where high concentrations of hydrogen 
sulfide, moisture, and oxygen are present in the atmosphere, the deterioration of concrete is 
caused mainly by biogenic acids. The so-called microbially-induced concrete corrosion in sewer 
systems has been a serious unsolved problcnQ for long time. The presence of microorganisms 
such as fungi, algae or bacteria can induce formation of aggressive biofilms on concrete surfaces. 
Particularly, the sulfuric acid that causes corrosion of sewer crowns is generated by such a complex 
microbial ecosystem especially in hot environments. The precise role of microorganisms in the 
context of sulfates attack on concrete (here we focus on sewer pipes) is quite complex and is 
therefore less understood from both experimental and theoretical points of view; see, e.g., the 
experimental studies [3] (optimum pH and growth kinetics of four relevant bacterial strains), 
1 11 j (characteristics of the crown microbial system), |12j (microbiologically influenced corrosion 
of natural sandstone), |16) (succession of sulfur-oxidizing bacteria in the bacterial community 
on corroding concrete), |17) (isolation of Thiobacillus thiooxidans) , [21] (Hamburg sewers) . (24] 
(air- water transfer of hydrogen sulfide). As a consequence of this, an accurate large-time forecast 
of the penetration of the sulfate corrosion front is very difficult to obtain. 
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1 There arc a lot of financial implications if you want to change the network of pipes in a city like Fukuoka. 
Our statement here is that questions like Why changing the pipes if corrosion is not so strong yet and therefore 
the mechanics structure of the network can/could still hold for 5 more years? can be addressed in a rigorous 
mathematical multiscale framework. Such an approach would allow a good understanding and prediction at least 
of extreme situations. 
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We want to stress the fact that concrete, in spite of its strong heterogeneity, is mechanically 
a well-understood material with known composition. Also, the cement (paste) chemistry is well 
understood. However, all cement-based materials (including concrete) involve a combination 
of "heterogeneous multi-phase material", "multiscale chemistry", "multiscale transport" (flow, 
diffusion, ionic fluxes, etc.), and "multiscale" mechanics. Having in view this complexity, such 
materials are sensu stricto very difficult to describe, to analyze mathematically, and last but not 
least, to deal with numerically. We expect that only after the multiscale aspects of such materials 
are handled properly, good predictions of the large-time behavior may be obtained. This is our 
path to addressing this corrosion scenario that is often referred to as the sulfatation problem. 

Before closing these background notes, let us add some remarks [5D] on a closely-related topic 
of acid sulphate soils^J which might attract the attention of the multiscale research in the near 
future. Acid sulphate soils are an important class of soils worldwide. Particularly in coastal 
areas, sediments often contain a large amount of iron sulfide (FeS and/or FeS2). When by 
drainage the sediment is exposed to air, this iron sulfide will oxidize to iron sulphate (FeSCU). 
As long as the sediment still contains calcium carbonate (G&CO^), the FeS will react with it, 
resulting in gypsum (CaSC>4) and iron oxide (Fe2 0a). Gypsum, being much more soluble than 
calcium carbonate, will tend to leach to the ground- and surface-waters. If the FesC>4 is no longer 
removed by reactions with CaC03 or other materials, it will tend to accumulate, resulting in a 
drop of the pH below 4. The problematic acid sulphate soil then will have become a reality. Acid 
sulphate soils were first described in 1886 by the Dutch Chemist Jacob Maarten van Bcmmelen, 
in connection with problems arising in the Haarlemmermeer Polder. Much of the work in this 
direction is/was done in the tropics, including Indonesia, Vietnam, and Australia; see, e.g., [13] . 

1.2. Objectives and structure of the paper. In order to be able to tackle the biophysics 
of the problem at a later stage, eventually coupled with the mechanics of the concrete and the 
actual capturing of the macroscopic fracture initiation, we focus here on a much simpler setting 
modeling the multiscale transport and reaction of the active chemical species involved in the 
sulfatation process. Therefore, the approach and results reported here are only preliminary. 

Our main objective is twofold: using a multiscale reaction-diffusion system for concrete cor- 
rosion (that allows for feedback between micro and macro scales), 

• calculate pH profiles and detect the eventual presence of "sudden" pH drops; 

• extract from gypsum concentration profiles the approximate position of macroscopic 
corrosion fronts. 

In Section[2j we present the reaction mechanisms taking place in sewer pipes. In Section[3j we 
give a mathematical description of the problem and we set a two-scale PDE-ODE system. We 
briefly comment on a few mathematical properties of the model. In Section [4j we approximate a 
macroscopic pH numerically using a multiscale FD scheme and comment briefly on the numerical 
results. 



2. A FEW NOTES ON THE INVOLVED CHEMISTRY 

Our model includes two important features: 

• continuous transfer of H2S from water to air phase and vice versa; 

• fast production of gypsum at solid-water interface. 

We incorporate the Henry's law to model the transfer of H2S from the water to the air phase 
and vice versa [2j 124j . The production of gypsum at the solid- water interface is modeled by a 
non-linear reaction rate, given by (15). 



2 Compared to concrete, soils are much easier to handle. Their mechanics is simpler and their chemistry is 
often rudimentary, if any. 
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Figure 1. Left: Cross-section of a sewer pipe. Middle: Mesoscopic periodic 
approximation of a REV. Right: Our concept of pore geometry (microstructure). 



There are many variants of severe attack to concrete in sewer pipes which influence the per- 
formance of concrete structure depending on the intensity of the reactions, the environment, and 
the turbulence of the wastewater |23j . We focus here on the most aggressive one, namely we 
consider the following reaction mechanisms causing sulfatation, viz. 

(1) 10 H + + SC>4~ + org. matter -> H 2 S(aq) + 4H 2 + oxidized matter 

(2) H 2 S(aq) + 2 2 -> 2H+ + SO^ 

(3) H 2 S(aq) ^ H 2 S(g) 

(4) 2 H 2 + H+ + SOl" + CaCC-3 -> CaS0 4 • 2 H 2 + HC0 3 _ 

Reaction ^ is typically a surface reaction taking place as soon as water and air phases meet 
together. It plays an important role in transferring the H 2 S from the air phase to the liquid phase 
where the corrosion actually takes place. For modeling details such as a Henry-like "reaction" 
mechanism, we refer the reader to |S] and references cited therein. 



3. MULTISCALE DESCRIPTION OF THE SULFATATION PROBLEM 

We assume that the geometry of our concrete sample (porous medium) consists of a system 
of pores periodically distributed inside a three-dimensional cube :— [a, b] 3 with a, b G E and 
b > a. The exterior boundary of 51 consists of two disjoint, sufficiently smooth parts: the 
Neumann boundary T N and the Dirichlet boundary T D . We assume that the pores in concrete 
are made of stationary water film, air and solid parts in different ratios depending on the local 
porosity. The reference pore, say Y := [0, l] 3 , has three pair-wise disjoint domains Y 3 , Y w and 
Y a with smooth boundaries T sw and T wa as shown in Fig. [T] such that 

Y = Y S L)Y W U Y a . 

We refer the reader to |S] for more description of the multiscale geometry of the porous 
material. For a single scale (macroscopic) approach of a sulfatation scenario, we refer the reader 
to pQ, e.g. 

We consider a two-scale system of PDEs and one ODE for unknown functions U\ : f2 x (0, T) — >• 
E, u k : fi x Y w x (0, T) -> E, k G {2, 3}, and u 4 : Q x T sw x (0, T) -> E where (0, T) is the time 
interval. The model under consideration is derived by formal homogenization using different 
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scalings of the diffusion coefficients in [5j (see also [H]) and is given by 

(5a) d t ux - efxAui = —B \ Hu x - J u 2 dj y ^j , in fi x (0, T), 

(5b) p 2 d t u 2 - p 2 d 2 A y u 2 = -<§> 2 k 2 u 2 + <S> 3 k 3 u 3 , infixTx (0, T), 

(5c) /3 3 a t M 3 - /3 3 d 3 A^3 = $ 2 fc 2 M2 - $ 3 fc 3 w 3 , infixTx (0, T), 

(5d) /3 4 <9 t u 4 = $ 4 M(w 3 ,w 4 ), in fi x r SUJ x (0,T), 

where ui denotes the concentration for H2S gaseous species, u 2 for H2S aqueous species, U3 
for H2SO4, and U4 for gypsum at r s, "\ The water film is taken here to be stationary. A 
detailed modeling of the role of water is still open, see, e.g., [IH1 131 [52]. A without subscript 
denotes the Laplace operator with respect to macroscopic variable x and A y with respect to 
microscopic variable y. d-f y represents the differential over the surface T wa . /3k > 0, k £ 
{2, 3, 4}, represents the ratio of the maximum concentration of the fc-th species to the maximum 
concentration of H2S0 4 , d. L > 0, i £ {1,2,3}, are the diffusion coefficients, B is a dimensionless 
Biot number which gives the mass transfer rate between water and air phases, and kj : Y — > 
i?, j £ {2,3,4}, are functions modeling the reaction rate "constants". <&k (k £ {2,3,4}) are 
Damkohler numbers corresponding to three distinct chemical mechanisms (reactions). They are 
dimensionless numbers comparing the characteristic time of the fastest transport mechanism 
(here, the diffusion of H2S in the gas phase) to the characteristic timescale of the fe-th chemical 
reaction. 

The system ^ is supplemented with initial and boundary conditions, which read as 

(6) ui(x,0) = ul(x), onf7x(0,T), 

(7) u k (x,y,0)=u° k (x,y), on fi x Y w x (0, T), k £ {2, 3}, 

(8) Ui {x, y, 0) - ul{x, y), on tlx T sw x (0, T), 

(9) ux=uf, onr D x(0,T), 

(10) n N ■ (diVui) = 0, on T N x (0, T), 

(11) n wa ■ (d 2 W y u 2 ) = B \Hui -J u 2 d 7 ^ , on £1 x T wa x (0, T), 

(12) n sw ■ (d 2 V y u 2 ) = 0, on fi x T sw x (0, T), 

(13) n wa ■ (d 3 V y u 3 ) = 0, on fi x T wa x (0, T), 

(14) n sw ■ (d 3 V y u 3 ) = -$ 3 ri(u3, Ui ), on fi x T sw x (0, T), 

where denotes the outward unit normal vector to 9fi along T N , and n wa and n sw denote the 
outward unit normal vectors to Y w along T wa and T sw , respectively. Note that the "information" 



at the micro-scale is connected to the macro-scale situation via the right-hand side of ( 5a ) and via 



the micro-macro boundary condition (111. It is also worth noticing that all involved parameters 
(except for H, d 3 and B) contain microscopic information. The coefficients d 3 and B are effective 
ones (see [9l [8] for a way of calculating them) , while H can be read off from existing macroscopic 
experimental data. 

We consider the following form of the reaction rate r\ at the interface T 3W 

(15) r,(a,0) 



a p (/3-/3)«, ifa>0,/3>0, 
0, otherwise, 



where (3 is a known maximum concentration of gypsum at T ew and p > l,q > 1 are partial 
orders of reaction. For more modeling possibilities of choosing 77, see [10] . It is worth noting that 
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production terms like 

are usually referred in the literature as Henry's or Raoult's law, where H > is known Henry's 
constant. 

We refer the reader to Theorem 3] for statements regarding the global existence and 
uniqueness a weak solution to problem (p^-(|T4|) (see also [15] for the analysis on a closely related 
problem). 

4. Simulation at a macroscopic pH scale. Capturing free boundaries 

In this section the model ([5| is applied to the simulation of the acid corrosion due to a 
microbiotical layer on a cement specimen. We focus on extracting the position of the corrosion 
front and on the acid reaction, which we use to obtain macro-scale profiles of pH. Both of these 
results can be compared to experimental data published, e.g., in 116j . 
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Figure 2. Time evolution of u\ (concentration of H2S(g)) shown at t € 
{2000, 4000, 8000, 12000, 16000, 20000} in left-right and top-bottom order. 
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Figure 3. Time evolution of it 4 (gypsum) shown at t £ 
{2000,4000,8000, 12000, 16000,20000} in left-right and top-bottom order. 

For the purpose of the simulation we employ a numerical scheme for a reduced 1D/2D version 
of the system The reduction consists in taking fl := (0, L) and Y w = (0, £) as one-dimensional 
intervals, which in effect corresponds to analysing the specimen in a perpendicular direction to its 
surface away from edges and to simplifying the micro cell geometry, respectively. The numerical 
scheme is based on the method of lines, where in space we use finite difference discretization 
and in time we employ an implicit higher-order time integrator for the solution of the non-linear 
ODE system. See [51 [7] for further details of the numerical scheme, its analysis with respect 
to convergence to the weak solutions and some basic numerical experiments. For details of a 
computer implementation of the numerical scheme we refer the reader to [141 Chapter 7] . 

In Table [l] we summarize values of the model parameters used in the simulations described 
below. 

4.1. Free boundaries. Figures [2] and [3] show the evolution of ui(x,t) and u±{x,t) in time. The 
Dirichlet boundary condition iii(0,i) = itf models a constant inflow of H2S(g) at x = 0, i.e., 
at the surface of the specimen. As the gas diffuses through the porous structure, it enters the 
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Figure 4. Position of the corrosion front. 
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Table 1. Parameter values used in the numerical simulation. 



water film in the pores due to the reaction ([3]) , where it undergoes biogenic oxidation to sulfuric 
acid. Consequently, its concentration decreases with increasing depth. As the system becomes 
saturated and as the sulfatation reaction Q converts available cement into gypsum, the total 
concentration of H 2 S(g) starts to increase (Fig. [2]). 

Sulfuric acid that arises from the oxidation of H2S(aq) then reacts at y = I with the cement 
paste and converts it into gypsum (u 4 ) whose concentration profile is shown in Fig. [3] Inter- 
estingly, although the behavior of u% is as expected (i.e., purely diffusive), we notice that a 
macroscopic gypsum layer (region where 1*4 is produced) is formed around t — 1500 and grows 
in time. The figure clearly indicates that there are two distinct regions separated by a slowly 
moving intermediate layer: the left region — the place where the gypsum production reached sat- 
uration (a threshold), and the right region — the place of the ongoing sulfatation reaction Q (the 
gypsum production has not yet reached here the natural threshold). 

We use U4 to extract an approximate position of the corrosion front p(t) which we define as 
(in our scenario, we expect U4 to be decreasing) 

p(t) :={x£ (0,L) I Ui{x,t) =$-e}, 

where s is a small parameter. Figure [4] shows a graph of p(t) arising from our numerical exper- 
iment. We notice that as the corroding front advances further into the concrete specimen, its 
rate of growth decreases. This is in agreement with experimental data since the hydrogen sulfide 
gas supplied from the outside environment has to be transported (by diffusion) over ever larger 
distance. It is important to note that the precise position of the separating layer is a priori 
unknown and to capture it simultaneously with the computation of the concentration profile 
would require a moving-boundary formulation similar to the one reported in [5]. 
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Figure 5. Time evolution of macro-scale pH profiles computed from micro-scale 
information shown at t € {2000, 4000, 8000, 12000, 16000, 20000} in left-right and 
top-bottom order. 



4.2. Drop in pH. Emission of hydrogen sulfide from the wastewater to the air space of sewer 
pipe is an important process because the problems of hydrogen sulfide in sewer pipes are associ- 
ated with gaseous hydrogen sulfide. Hydrogen sulfide is a weak acid with a dissociation constant 
of 7.0 (at 20°C) and only the non-dissociated form is emitted in the air space sewer pipe. The 
pH of the wastewater is therefore of importance when evaluating the potential hydrogen sulfide 
emission. After the hydrogen sulfide arrives at and diffuses into the concrete, the oxidation of 
hydrogen sulfide is biological once the pH of the solid matrix has dropped below approximately 
8-9 [18] . This represents the tendency of hydronium ions to interact with other components of 
the solution, which affects among other things the electrical potential read using a pH meter. 
The concentration of hydrogen ions is expressed as pH scale and pH is defined as a negative 
decimal logarithm of the concentration of hydronium ions dissolved in a solution. 
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As the model considered in this paper contains information both from micro (pore) scale and 
macro-scale, we need a way of computing macro-scale values of pH from the available micro- 
scale data. Such information is not readily available in the system ([5]). However, sulfuric acid 
is a diprotic acid with two stages of dissociation, where the first stage occurs fully and the 
dissociation in the second stage can be neglected. Therefore, the concentration of hydronium 
ions is proportional to the concentration of sulfuric acid, which is represented by 113 in our model. 
We extract the macro-scale concentration of sulfuric acid at each x by taking a volume average 
of M3 over Y w . Thus, we use the following expression for computing macroscopic pH: 

(!6) pH macro (x, t) = - log 10 f T^jrf J y u s(x, V, t) dy 

where k a is the activity of hydronium ions. 

The macro-scale pH profile computed using the formula ( 16 1 is shown in Figure U We can see 



that at the beginning of the simulation (first graph) with increasing depth the pH also increases 
from acid to more basic values as expected. Once all the available cement is consumed and 
converted into gypsum (this happens for the first time at x = between the first and second 
graph in Figure [5] around t = 1500), the pH drops rapidly across the corrosion front. This is 
due to the fact that behind the corrosion front the sulfuric acid is no longer neutralised by the 
sulfatation reaction Q. 

Note that our pH profiles are not in the experimental range. We expect the size of the 
drop will become comparable to the one seen in experiments as soon as effects of nonlinear 
moisture transport and bacteria motility and chemical activity are taken into account in the 
model equations. The main message that we want to draw is that we are able to detect and 
compute a macroscopic pH drop, once the needed micro-information is available. 
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